Travis continuation integration Test Cases

In [1]:
import numpy as np
import os
import pandas as pd
import itertools
from PIL import Image
import seaborn as sns
import fitsio
import skimage.io
import galsim

from astrometry.util.fits import fits_table, merge_tables


# to make this notebook's output stable across runs
np.random.seed(7)

# To plot pretty figures
#%matplotlib inline
%matplotlib notebook

import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['axes.labelsize'] = 14
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12

%load_ext autoreload
%autoreload 2
In [2]:
from obiwan.qa.visual import readImage,sliceImage,plotImage, flux2mag
In [3]:
# Testcase functions not in py/obiwan/ directory
a=os.path.join(os.getcwd(),
              "../../tests/end_to_end")
!ls $a

import sys
sys.path.append(a)
__init__.py                             out_testcase_decals_z_dr5_elg_onedge
legacypipedir_1238p245_dataset_DR3      out_testcase_decals_z_dr5_star
legacypipedir_1238p245_dataset_DR5      outdir_1238p245_DR3
out_testcase_bassmzls_grz_dr6_elg       outdir_1238p245_DR5
out_testcase_decals_grz_dr5_elg         test_datasets.py
out_testcase_decals_z_dr5_elg           testcase_bassmzls_grz
out_testcase_decals_z_dr5_elg_allblobs  testcase_decals_grz
out_testcase_decals_z_dr5_elg_coadds    testcase_decals_z
In [5]:
print('hello from docker image')
hello from docker image

testcase_DR_z (z-band, 200x200 pix region, real galaxy at center)

all_blobs == True (fit model to all detected sources)

In [17]:
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=False,all_blobs=True)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_elg_allblobs/elg
In [13]:
fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(t.img_fits['z'],ax,qs=None)
plotImage().circles(t.simcat.x,t.simcat.y,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=5./0.262,color='y')
plotImage().circles(t.obitractor[ireal].bx,t.obitractor[ireal].by,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=4./0.262,color='m')
In [18]:
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y),
                                   r=3.5/t.brickwcs.pixel_scale())
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])
apy_table = photutils.aperture_photometry(t.sims_fits['z'], apers)
sims_apflux= np.array(apy_table['aperture_sum'])

sky_apflux= img_apflux - sims_apflux
#sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

print('simcat id',t.simcat.id)
print('sky_apflux:',sky_apflux)
#print('sky_meas:',sky_meas)
print('sims_apflux:',sims_apflux)
print('simcat_flux:',t.simcat.zflux)
print('obitrac_flux:',t.obitractor[itrac].flux_z)
print('img_apflux:',img_apflux)
print('obitrac_apflux:',t.obitractor[itrac].apflux_z[:,5])

print(sky_apflux + sims_apflux - t.obitractor[itrac].apflux_z[:,5])
print(t.obitractor[itrac].flux_z - t.simcat.zflux)
('simcat id', array([1, 2, 3, 4]))
('sky_apflux:', array([-0.8440134 , -0.56180657,  0.44062445, -1.00573623]))
('sims_apflux:', array([ 36.60767405,  36.59797994,  36.61960024,  36.56948628]))
('simcat_flux:', array([ 38.28181781,  38.28181781,  38.28181781,  38.28181781]))
('obitrac_flux:', array([ 33.39308548,  32.36734772,  34.74517059,  32.90026855], dtype=float32))
('img_apflux:', array([ 35.76366065,  36.03617336,  37.06022469,  35.56375005]))
('obitrac_apflux:', array([ 35.9729805 ,  36.11907196,  36.9563179 ,  35.5991478 ], dtype=float32))
[-0.20931985 -0.0828986   0.10390679 -0.03539775]
[-4.88873233 -5.91447009 -3.53664721 -5.38154925]
In [24]:
np.max([t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapeexp_r],axis=1)
Out[24]:
array([ 0.44706032,  0.44706032], dtype=float32)
In [19]:
t.obitractor[ireal].apflux_z[:,5], t.obitractor[ireal].flux_z
Out[19]:
(array([ 38.94617462], dtype=float32), array([ 38.28445053], dtype=float32))
In [20]:
t.simcat.rhalf,t.simcat.e1,t.simcat.e2,t.simcat.n
Out[20]:
(array([ 0.37005,  0.37005,  0.37005,  0.37005]),
 array([ 0.33340144,  0.33340144,  0.33340144,  0.33340144]),
 array([ 0.27338931,  0.27338931,  0.27338931,  0.27338931]),
 array([4, 4, 4, 4]))
In [21]:
t.obitractor[itrac].type,t.obitractor[itrac].shapeexp_r
Out[21]:
(array(['REX ', 'REX ', 'REX ', 'REX '],
       dtype='|S4'),
 array([ 0.44169348,  0.38656846,  0.44706032,  0.40121055], dtype=float32))
In [22]:
t.obitractor[itrac].shapeexp_r-t.simcat.rhalf
Out[22]:
array([ 0.07164348,  0.01651846,  0.07701032,  0.03116055])
In [11]:
t.obitractor[itrac].shapeexp_r-t.simcat.rhalf
Out[11]:
array([ 0.0719691 ,  0.01574509,  0.07740407,  0.03068192])
In [6]:
def flux2mag(flux):
    return -2.5*np.log10(1e-9 * flux)

def mag2flux(mag):
    return 10**(9-mag/2.5)

t.obitractor[ireal].type,t.obitractor[ireal].flux_z, flux2mag(t.obitractor[ireal].flux_z)
Out[6]:
(array(['DEV '],
       dtype='|S4'),
 array([ 38.28407288], dtype=float32),
 array([ 18.54245377], dtype=float32))
In [50]:
t.obitractor[ireal].shapedev_r,t.obitractor[ireal].shapedev_e1,t.obitractor[ireal].shapedev_e2
Out[50]:
(array([ 0.37025318], dtype=float32),
 array([ 0.33367831], dtype=float32),
 array([ 0.27334914], dtype=float32))
In [51]:
print(t.obitractor[itrac].type,t.obitractor[itrac].flux_z, flux2mag(t.obitractor[itrac].flux_z))
t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapeexp_e1,t.obitractor[itrac].shapeexp_e2
(array(['REX ', 'REX ', 'REX ', 'REX '],
      dtype='|S4'), array([ 33.38793564,  32.36148453,  34.74031448,  32.89451981], dtype=float32), array([ 18.69102669,  18.72492981,  18.64791489,  18.70719147], dtype=float32))
Out[51]:
(array([ 0.44197041,  0.38681737,  0.44734058,  0.40145969], dtype=float32),
 array([ 0.,  0.,  0.,  0.], dtype=float32),
 array([ 0.,  0.,  0.,  0.], dtype=float32))
In [31]:
x=np.array([198,2,100,150]) + t.zoom[0]
y=np.array([2,198,2,198]) + t.zoom[2]
t.brickwcs.pixelxy2radec(x,y)
Out[31]:
(array([ 174.23933751,  174.25500639,  174.24716429,  174.24318502]),
 array([ 24.32087412,  24.33512531,  24.32086773,  24.33513542]))
In [47]:
t.simcat.ra
Out[47]:
array([ 174.24714385,  174.24715649,  174.25116781,  174.24313254])

all_blobs (coadds, top to bottom: image, simulated image, and residual)

In [49]:
fig,ax=plt.subplots(3,3,figsize=(7,6))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])

    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.obitractor.bx,t.obitractor.by,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=5./0.262,color='y')
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')

    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])
In [12]:
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(t.blobs,ax[0,0],qs=None)
plotImage().imshow(t.img_jpg,ax[0,1],qs=None)
plotImage().imshow(t.model_jpg,ax[1,0],qs=None)
plotImage().imshow(t.resid_jpg,ax[1,1],qs=None)
In [15]:
plot3d(t.ivar_fits['z'])
In [50]:
len(isim),len(itrac),len(ireal)
Out[50]:
(4, 4, 1)
In [51]:
t.simcat.zflux, t.obitractor[itrac].flux_z
Out[51]:
(array([ 38.26484866,  38.26484866,  38.26484866,  38.26484866]),
 array([ 35.25556183,  34.26049423,  36.58589172,  34.72808075], dtype=float32))
In [52]:
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y),
                                   r=3.5/t.brickwcs.pixel_scale())
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])
apy_table = photutils.aperture_photometry(t.sims_fits['z'], apers)
sims_apflux= np.array(apy_table['aperture_sum'])

sky_apflux= img_apflux - sims_apflux
#sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

print('sky_apflux:',sky_apflux)
#print('sky_meas:',sky_meas)
print('sims_apflux:',sims_apflux)
print('simcat_flux:',t.simcat.zflux)
print('obitrac_flux:',t.obitractor[itrac].flux_z)
print('obitrac_apflux:',t.obitractor[itrac].apflux_z[:,5])


# apers= photutils.CircularAperture(([100],[100]),
#                                    r=3.5/brickwcs.pixel_scale())
# apy_table = photutils.aperture_photometry(img_fits['g'], apers)
# real_apflux= np.array(apy_table['aperture_sum'])

# print(simcat.gflux,
#       sim_apflux,
#       real_apflux)
# print(flux2mag(simcat.gflux) - flux2mag(sim_apflux),flux2mag(simcat.gflux[0]),flux2mag(real_apflux))

('sky_apflux:', array([-0.84267303, -0.56008409,  0.44084312, -1.00828833]))
('sims_apflux:', array([ 38.55891982,  38.5594253 ,  38.56627197,  38.53363932]))
('simcat_flux:', array([ 38.26484866,  38.26484866,  38.26484866,  38.26484866]))
('obitrac_flux:', array([ 35.25556183,  34.26049423,  36.58589172,  34.72808075], dtype=float32))
('obitrac_apflux:', array([ 37.90312195,  38.06328583,  38.90242767,  37.55826569], dtype=float32))
In [53]:
print(sky_apflux + sims_apflux - t.obitractor[itrac].apflux_z[:,5])
[-0.18687516 -0.06394461  0.10468742 -0.03291469]
In [54]:
print(t.simcat.zflux + sky_apflux - t.obitractor[itrac].apflux_z[:,5])
[-0.48094632 -0.35852125 -0.19673589 -0.30170536]
In [55]:
apers= photutils.CircularAperture((t.obitractor[ireal].bx,t.obitractor[ireal].by),
                                   r=3.5/t.brickwcs.pixel_scale())
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])

print('real galaxy apflux:',img_apflux)
print('real galaxy obitrac flux:',t.obitractor[ireal].flux_z)
print('real galaxy obitrac apflux:',t.obitractor[ireal].apflux_z[0][5])
print(flux2mag(img_apflux) - flux2mag(t.obitractor[ireal].flux_z))
('real galaxy apflux:', array([ 38.91102173]))
('real galaxy obitrac flux:', array([ 38.26608276], dtype=float32))
('real galaxy obitrac apflux:', 38.91172)
[-0.01814652]
In [9]:
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=True)
import plotly.graph_objs as go

def plot3d(image,zlim=None):
    data = [
        go.Surface(
            z=image
        )
    ]
    layout=go.Layout()
    if zlim:
        layout = go.Layout(
                        scene= dict(
                            zaxis = dict(
                                nticks=10, range = [zlim[0],zlim[1]]),
                        )
                      )
    fig = go.Figure(data=data, layout=layout)
    #fig = go.Figure(data=data)
    iplot(fig, filename='image')


In [ ]:
plot3d(img_fits)
In [ ]:
plot3d(sims_fits)
In [24]:
plot3d(img_fits - sims_fits)

Early coadds (just make the coadds, don’t fit models)

In [4]:
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_grz',dataset='DR5',
             all_blobs=True,early_coadds=True)

t= AnalyzeTestcase(**kwargs)
t.load_outputs()
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_grz_elg_allblobs_coadds/elg/028/0285m165/rs0/
In [5]:
fig,ax=plt.subplots(3,3,figsize=(7,6))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])

    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')

    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])

stars (inject stars into the images)

In [60]:
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             obj='star',
             onedge=False,add_noise=False,all_blobs=True)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_star_allblobs/star/174/1741p242/rs0/
In [61]:
fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(t.img_fits['z'],ax,qs=None)
plotImage().circles(t.simcat.x,t.simcat.y,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=5./0.262,color='y')
plotImage().circles(t.obitractor[ireal].bx,t.obitractor[ireal].by,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=4./0.262,color='m')
In [62]:
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y),
                                   r=3.5/t.brickwcs.pixel_scale())
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])
apy_table = photutils.aperture_photometry(t.sims_fits['z'], apers)
sims_apflux= np.array(apy_table['aperture_sum'])

sky_apflux= img_apflux - sims_apflux
#sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

print(t.simcat.id)
print('sky_apflux:',sky_apflux)
#print('sky_meas:',sky_meas)
print('sims_apflux:',sims_apflux)
print('simcat_flux:',t.simcat.zflux)
print('obitrac_flux:',t.obitractor[itrac].flux_z)
print('obitrac_apflux:',t.obitractor[itrac].apflux_z[:,5])
added_flux= np.array([39.135,39.102,39.122,39.115])
print('added_flux/1.02',added_flux/1.02)
[1 2 3 4]
('sky_apflux:', array([-0.84315645, -0.56062535,  0.4417466 , -1.00764566]))
('sims_apflux:', array([ 37.83162486,  37.85237236,  37.85502317,  37.80914618]))
('simcat_flux:', array([ 38.26484866,  38.26484866,  38.26484866,  38.26484866]))
('obitrac_flux:', array([ 36.89432526,  37.22283554,  37.80667496,  37.30305481], dtype=float32))
('obitrac_apflux:', array([ 37.17383194,  37.24513245,  38.22380829,  36.844944  ], dtype=float32))
('added_flux/1.02', array([ 38.36764706,  38.33529412,  38.35490196,  38.34803922]))
In [8]:
t.obitractor[itrac].type
Out[8]:
array(['PSF ', 'PSF ', 'PSF ', 'PSF '],
      dtype='|S4')
In [12]:
dflux= t.simcat.zflux + \
        sky_apflux - t.obitractor[itrac].apflux_z[:,5]
print('diff flux: ',dflux)
dmag= flux2mag(t.simcat.zflux + sky_apflux) - flux2mag(t.obitractor[itrac].apflux_z[:,5])
print('diff mag: ',dmag)
('diff flux: ', array([-0.61297537, -0.36756682, -0.36557727, -0.42651694]))
('diff mag: ', array([ 0.01764039,  0.01053395,  0.01020564,  0.01235859]))

onedge (inject galaxies on the edges of the images, how good is best-fit model?)

In [19]:
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             onedge=True,add_noise=False,all_blobs=False)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()

print(t.simcat.id[isim])
print([(x,y) for x,y in zip(t.simcat.x[isim],t.simcat.y[isim])])
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_elg_onedge/elg
rhalf [-0.01998445  0.13738838  0.02882487  0.11826912]
apflux [ 0.1496385   0.09753577 -0.00312189  0.01799367]
skyflux [-0.26486025 -0.06903909  1.73439526 -0.05648272]
modelflux [ 0.33289523 -0.195715    5.08362384  2.63030047]
passed
[1 2 3 4]
[(197.99994482662373, 2.0000151657400238), (1.9999776837685204, 198.00002486870926), (100.00003559541528, 2.0000627804956821), (150.00001340770677, 197.9999763324513)]
In [20]:
fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(t.img_fits['z'],ax,qs=None)
plotImage().circles(t.simcat.x,t.simcat.y,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=5./0.262,color='y')
plotImage().circles(t.obitractor[itrac].bx,t.obitractor[itrac].by,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=4./0.262,color='m')


In [18]:
rhalf= np.max((t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapedev_r),
              axis=0)
e1= np.max((t.obitractor[itrac].shapeexp_e1,t.obitractor[itrac].shapedev_e1),
              axis=0)
e2= np.max((t.obitractor[itrac].shapeexp_e2,t.obitractor[itrac].shapedev_e2),
              axis=0)
print(t.simcat[isim].id)
print(t.obitractor[itrac].type,rhalf,e1,e2)
print(t.simcat.n[isim],t.simcat.rhalf[isim],t.simcat.e1[isim],t.simcat.e2[isim])
print(t.simcat.zflux[isim] - t.obitractor.flux_z[itrac])
[1 2 3 4]
(array(['DEV ', 'REX ', 'REX ', 'REX '],
      dtype='|S4'), array([ 0.34959555,  0.50696838,  0.39840487,  0.48784912], dtype=float32), array([ 0.11206593,  0.        ,  0.        ,  0.        ], dtype=float32), array([ 0.0281052,  0.       ,  0.       ,  0.       ], dtype=float32))
(array([4, 4, 4, 4]), array([ 0.36958,  0.36958,  0.36958,  0.36958]), array([ 0.33375,  0.33375,  0.33375,  0.33375]), array([ 0.27345,  0.27345,  0.27345,  0.27345]))
[ 0.33289523 -0.195715    5.08362384  2.63030047]
In [24]:
import galsim
a=galsim.PositionD(10,10)
a.x

Out[24]:
10.0
In [39]:
a=galsim.Image(np.zeros((10,10)))
a.shift(dx=0.5,dy=0.5)
a.setCenter(10,10)
In [40]:
a.trueCenter(),a.center()
Out[40]:
(galsim.PositionD(x=9.5, y=9.5), galsim.PositionI(x=10, y=10))

all_blobs == False (only fit models inside blobs that have fake sources)

In [22]:
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=False,all_blobs=False)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_elg/elg

The central galaxy is real and was not modeled

In [25]:
fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(t.img_fits['z'],ax,qs=None)
plotImage().circles(t.simcat.x,t.simcat.y,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=7./0.25,color='y')
plotImage().circles(t.obitractor[itrac].bx,t.obitractor[itrac].by,ax,
                    img_shape=t.img_fits['z'].shape,
                    r_pixels=4./0.25,color='m')
In [23]:
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y),
                                   r=3.5/t.brickwcs.pixel_scale())
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])
apy_table = photutils.aperture_photometry(t.sims_fits['z'], apers)
sims_apflux= np.array(apy_table['aperture_sum'])

sky_apflux= img_apflux - sims_apflux
#sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

print('sky_apflux:',sky_apflux)
#print('sky_meas:',sky_meas)
print('sims_apflux:',sims_apflux)
print('simcat_flux:',t.simcat.zflux)
print('obitrac_flux:',t.obitractor[itrac].flux_z)
print('obitrac_apflux:',t.obitractor[itrac].apflux_z[:,5])
('sky_apflux:', array([-0.8424493 , -0.56035691,  0.44068548, -1.00731097]))
('sims_apflux:', array([ 35.84701954,  35.87222211,  35.87172834,  35.8296502 ]))
('simcat_flux:', array([ 38.28181781,  38.28181781,  38.28181781,  38.28181781]))
('obitrac_flux:', array([ 32.69493866,  31.70311165,  34.05732346,  32.2202034 ], dtype=float32))
('obitrac_apflux:', array([ 35.21274948,  35.39452744,  36.20796204,  34.85829544], dtype=float32))
In [24]:
print(t.simcat.rhalf[isim],t.simcat.e1[isim],t.simcat.e2[isim])
print(t.obitractor.type[itrac],t.obitractor.shapeexp_r[itrac],t.obitractor.shapeexp_e1[itrac],t.obitractor.shapeexp_e2[itrac])
(array([ 0.37005,  0.37005,  0.37005,  0.37005]), array([ 0.33340144,  0.33340144,  0.33340144,  0.33340144]), array([ 0.27338931,  0.27338931,  0.27338931,  0.27338931]))
(array(['REX ', 'REX ', 'REX ', 'REX '],
      dtype='|S4'), array([ 0.4420191 ,  0.38579509,  0.44745407,  0.40073192], dtype=float32), array([ 0.,  0.,  0.,  0.], dtype=float32), array([ 0.,  0.,  0.,  0.], dtype=float32))
In [26]:
fig,ax=plt.subplots(3,3,figsize=(7,6))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])

    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.obitractor.bx,t.obitractor.by,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=5./0.262,color='y')
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')

    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])
In [11]:
plot3d(t.img_fits['z'])
In [ ]:
all_blobs == False
In [23]:
### No noise
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=False,all_blobs=False)
no= AnalyzeTestcase(**kwargs)
no.load_outputs()
isim,itrac,ireal= no.match_simcat_tractor()

print(no.obitractor[itrac].apflux_z)[:,5]
plot3d(no.sims_fits['z'])
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z/elg/174/1741p242/rs0/
[ 37.90312195  38.06328583  38.90242767  37.55826569]
In [24]:
### w/noise original formula
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=True,all_blobs=False)
noise= AnalyzeTestcase(**kwargs)
noise.load_outputs()
isim,itrac,ireal= noise.match_simcat_tractor()

print(noise.obitractor[itrac].apflux_z)[:,5]
plot3d(noise.sims_fits['z'])
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_addnoise/elg/174/1741p242/rs0/
[ 37.78057098  38.04665375  38.87749863  37.56063843]
In [40]:
### w/noise mine
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_z',dataset='DR5',
             add_noise=True,all_blobs=False)
my= AnalyzeTestcase(**kwargs)
my.load_outputs()
isim,itrac,ireal= my.match_simcat_tractor()

print(my.obitractor[itrac].apflux_z)[:,5]
plot3d(my.img_fits['z'])
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_z_addnoise/elg/174/1741p242/rs0/
[ 37.59292221  37.92685318  38.91709137  37.45846176]
In [31]:
sns.distplot(np.random.randn(100,100).flatten())
Out[31]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fce62556150>
In [8]:
print(no.obitractor[itrac].apflux_z[:,5] - my.obitractor[itrac].apflux_z[:,5])
print(no.obitractor[itrac].apflux_z[:,5] - noise.obitractor[itrac].apflux_z[:,5])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-566c0d0c5800> in <module>()
----> 1 print(no.obitractor[itrac].apflux_z[:,5] - my.obitractor[itrac].apflux_z[:,5])
      2 print(no.obitractor[itrac].apflux_z[:,5] - noise.obitractor[itrac].apflux_z[:,5])

NameError: name 'no' is not defined

5 galaxies but only 4 blobs, meaning only the 4 fake ones were modeled

In [62]:
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(t.blobs,ax[0,0],qs=None)
plotImage().imshow(t.img_jpg,ax[0,1],qs=None)
plotImage().imshow(t.model_jpg,ax[1,0],qs=None)
plotImage().imshow(t.resid_jpg,ax[1,1],qs=None)

testcase_DR_grz (Same as previous testcase, but three bands (grz) and a different region)

allblobs == True

In [83]:
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_grz',dataset='DR5',
             obj='elg',
             add_noise=False,all_blobs=True)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_grz_elg_allblobs/elg/028/0285m165/rs0/

Many sources are modeled, not just the 4 fake galaxies near center

In [84]:
fig,ax=plt.subplots(4,3,figsize=(7,8))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])

    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.obitractor.bx,t.obitractor.by,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=5./0.262,color='y')
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')

    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])
    plotImage().imshow(t.ivar_fits[b],ax[3,i])
In [23]:
plot3d(t.ivar_fits['r'])
In [85]:
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(t.blobs,ax[0,0],qs=None)
plotImage().imshow(t.img_jpg,ax[0,1],qs=None)
plotImage().imshow(t.model_jpg,ax[1,0],qs=None)
plotImage().imshow(t.resid_jpg,ax[1,1],qs=None)
In [71]:
len(isim),len(itrac),len(ireal)
Out[71]:
(4, 4, 1)
In [72]:
t.simcat.zflux, t.obitractor[itrac].flux_z
Out[72]:
(array([ 52.83479175,  52.83479175,  52.83479175,  52.83479175]),
 array([ 76.45344543,  74.84441376,  78.35695648,  77.82241821], dtype=float32))
In [78]:
import photutils
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y),
                                   r=3.5/t.brickwcs.pixel_scale())
for band in t.bands:
    print('band=',band)
    apy_table = photutils.aperture_photometry(t.img_fits[b], apers)
    img_apflux= np.array(apy_table['aperture_sum'])
    apy_table = photutils.aperture_photometry(t.sims_fits[b], apers)
    sims_apflux= np.array(apy_table['aperture_sum'])

    sky_apflux= img_apflux - sims_apflux
    #sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

    simcat_flux= t.simcat.get(b+'flux')
    trac_flux= t.obitractor[itrac].get('flux_'+b)
    trac_apflux=t.obitractor[itrac].get('apflux_'+b)[:,5]

    print('sky_apflux:',sky_apflux)
    #print('sky_meas:',sky_meas)
    print('sims_apflux:',sims_apflux)
    print('simcat_flux:',simcat_flux)
    print('obitrac_flux:',trac_flux)
    print('obitrac_apflux:',trac_apflux)
    print(sky_apflux + sims_apflux - trac_apflux)
    print(simcat_flux + sky_apflux - trac_apflux)
('band=', 'g')
('sky_apflux:', array([ 0.11542435,  0.0353508 ,  1.73446525,  0.74227955]))
('sims_apflux:', array([ 41.97473241,  42.02550686,  41.96238194,  42.05413268]))
('simcat_flux:', array([ 56.64281598,  56.64281598,  56.64281598,  56.64281598]))
('obitrac_flux:', array([ 52.03792953,  51.14485168,  54.86248016,  52.44288254], dtype=float32))
('obitrac_apflux:', array([ 41.97197723,  41.95433807,  43.69890976,  42.63935089], dtype=float32))
[ 0.11817952  0.10651959 -0.00206257  0.15706134]
[ 14.78626309  14.72382871  14.67837146  14.74574464]
('band=', 'r')
('sky_apflux:', array([ 0.11542435,  0.0353508 ,  1.73446525,  0.74227955]))
('sims_apflux:', array([ 41.97473241,  42.02550686,  41.96238194,  42.05413268]))
('simcat_flux:', array([ 56.64281598,  56.64281598,  56.64281598,  56.64281598]))
('obitrac_flux:', array([ 52.03792953,  51.14485168,  54.86248016,  52.44288254], dtype=float32))
('obitrac_apflux:', array([ 41.97197723,  41.95433807,  43.69890976,  42.63935089], dtype=float32))
[ 0.11817952  0.10651959 -0.00206257  0.15706134]
[ 14.78626309  14.72382871  14.67837146  14.74574464]
('band=', 'z')
('sky_apflux:', array([ 0.11542435,  0.0353508 ,  1.73446525,  0.74227955]))
('sims_apflux:', array([ 41.97473241,  42.02550686,  41.96238194,  42.05413268]))
('simcat_flux:', array([ 56.64281598,  56.64281598,  56.64281598,  56.64281598]))
('obitrac_flux:', array([ 52.03792953,  51.14485168,  54.86248016,  52.44288254], dtype=float32))
('obitrac_apflux:', array([ 41.97197723,  41.95433807,  43.69890976,  42.63935089], dtype=float32))
[ 0.11817952  0.10651959 -0.00206257  0.15706134]
[ 14.78626309  14.72382871  14.67837146  14.74574464]
In [74]:
t.img_fits['g'].shape
Out[74]:
(200, 200)
In [79]:
t.obitractor[ireal].type,t.obitractor[ireal].shapeexp_r,t.obitractor[ireal].shapedev_r,
Out[79]:
(array(['COMP'],
       dtype='|S4'),
 array([ 2.76125026], dtype=float32),
 array([ 0.24888119], dtype=float32))
In [80]:
for mod in ['shapedev','shapeexp']:
    print(t.obitractor[ireal].get(mod+'_e1'),t.obitractor[ireal].get(mod+'_e2'))
(array([ 0.22466888], dtype=float32), array([-0.1372889], dtype=float32))
(array([ 0.12009902], dtype=float32), array([ 0.2517896], dtype=float32))
In [81]:
for b in 'grz':
    print(b,flux2mag(t.obitractor[ireal].get('flux_'+b)))
('g', array([ 19.44617081], dtype=float32))
('r', array([ 18.69445229], dtype=float32))
('z', array([ 18.13280869], dtype=float32))
In [86]:
t.obitractor[itrac].type,t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapedev_r,
Out[86]:
(array(['EXP ', 'EXP ', 'DEV ', 'DEV '],
       dtype='|S4'),
 array([ 2.30459976,  2.32127619,  0.        ,  0.        ], dtype=float32),
 array([ 0.        ,  0.        ,  2.61173296,  2.21803832], dtype=float32))
In [87]:
t.obitractor[itrac].shapeexp_e1,t.obitractor[itrac].shapeexp_e2
Out[87]:
(array([ 0.06568176,  0.06891665,  0.        ,  0.        ], dtype=float32),
 array([-0.13783246, -0.14058098,  0.        ,  0.        ], dtype=float32))
In [88]:
t.obitractor[itrac].shapedev_e1,t.obitractor[itrac].shapedev_e2
Out[88]:
(array([ 0.        ,  0.        ,  0.11277404,  0.05932122], dtype=float32),
 array([ 0.        ,  0.        , -0.12263217, -0.13615023], dtype=float32))
In [90]:
for b in 'grz':
    print(b,flux2mag(t.obitractor[itrac].get('flux_'+b)))
('g', array([ 19.4939785 ,  19.47326469,  19.43455124,  19.55050659], dtype=float32))
('r', array([ 18.72727394,  18.7271843 ,  18.68118286,  18.78072548], dtype=float32))
('z', array([ 18.16558075,  18.18332863,  18.09988785,  18.22354126], dtype=float32))
In [66]:
apers= photutils.CircularAperture((t.obitractor[ireal].bx,t.obitractor[ireal].by),
                                   r=3.5/t.brickwcs.pixel_scale())
apy_table = photutils.aperture_photometry(t.img_fits['z'], apers)
img_apflux= np.array(apy_table['aperture_sum'])

print('real galaxy apflux:',img_apflux)
print('real galaxy obitrac flux:',t.obitractor[ireal].flux_z)
print('real galaxy obitrac apflux:',t.obitractor[ireal].apflux_z[0][5])
('real galaxy apflux:', array([ 40.51106904]))
('real galaxy obitrac flux:', array([ 56.64278793], dtype=float32))
('real galaxy obitrac apflux:', 40.513737)

allblobs == False

In [41]:
from test_datasets import Testcase,AnalyzeTestcase
kwargs= dict(name='testcase_DR5_grz',dataset='DR5',
             obj='elg',
             add_noise=False,all_blobs=False)
t= AnalyzeTestcase(**kwargs)
t.load_outputs()
t.simcat_xy()
#t.numeric_tests()
isim,itrac,ireal= t.match_simcat_tractor()
Loading from /home/kaylan/myrepo/obiwan/tests/end_to_end/out_testcase_DR5_grz_elg/elg

blobs are more extended, so less sources modeled, but still more than the 4 fake galaxies

In [42]:
fig,ax=plt.subplots(4,3,figsize=(7,8))
for i,b in enumerate(t.bands) :
    plotImage().imshow(t.img_fits[b],ax[0,i])

    plotImage().imshow(t.sims_fits[b],ax[1,i],qs=None)
    plotImage().circles(t.obitractor.bx,t.obitractor.by,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=5./0.262,color='y')
    plotImage().circles(t.simcat.x,t.simcat.y,ax[1,i],
                        img_shape=t.sims_fits[b].shape,
                        r_pixels=4./0.262,color='m')

    plotImage().imshow(t.img_fits[b] - t.sims_fits[b],ax[2,i])
    plotImage().imshow(t.ivar_fits[b],ax[3,i])
In [44]:
print(t.obitractor[itrac].type,t.obitractor[itrac].shapeexp_r,t.obitractor[itrac].shapedev_r)
print(t.simcat[isim].rhalf)
(array(['EXP ', 'EXP ', 'DEV ', 'DEV '],
      dtype='|S4'), array([ 2.30231023,  2.31898761,  0.        ,  0.        ], dtype=float32), array([ 0.        ,  0.        ,  2.60479164,  2.21619177], dtype=float32))
[ 2.  2.  2.  2.]
In [45]:
print(t.obitractor[itrac].shapeexp_e1,t.obitractor[itrac].shapedev_e1)
print(t.simcat[isim].e1)
(array([ 0.06559182,  0.06887048,  0.        ,  0.        ], dtype=float32), array([ 0.        ,  0.        ,  0.11189053,  0.05938281], dtype=float32))
[ 0.12408947  0.12408947  0.12408947  0.12408947]
In [46]:
print(t.obitractor[itrac].shapeexp_e2,t.obitractor[itrac].shapedev_e2)
print(t.simcat[isim].e2)
(array([-0.13772453, -0.14047767,  0.        ,  0.        ], dtype=float32), array([ 0.        ,  0.        , -0.12285928, -0.13606296], dtype=float32))
[ 0.26107663  0.26107663  0.26107663  0.26107663]
In [104]:
for b in 'grz':
    print(b,flux2mag(t.obitractor[itrac].get('flux_'+b)))
('g', array([ 19.49403381,  19.47325134,  19.43470192,  19.55048943], dtype=float32))
('r', array([ 18.72727203,  18.72712708,  18.68135262,  18.78062248], dtype=float32))
('z', array([ 18.16565132,  18.183321  ,  18.10007477,  18.22348976], dtype=float32))
In [28]:
plot3d(t.ivar_fits['g'])

only blobs containing fake sources are modeled, but this is still many sources

In [115]:
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(t.blobs,ax[0,0],qs=None)
plotImage().imshow(t.img_jpg,ax[0,1],qs=None)
plotImage().imshow(t.model_jpg,ax[1,0],qs=None)
plotImage().imshow(t.resid_jpg,ax[1,1],qs=None)
In [90]:
apers= photutils.CircularAperture((t.simcat.x,t.simcat.y),
                                   r=3.5/t.brickwcs.pixel_scale())
for band in t.bands:
    print('band=',band)
    apy_table = photutils.aperture_photometry(t.img_fits[b], apers)
    img_apflux= np.array(apy_table['aperture_sum'])
    apy_table = photutils.aperture_photometry(t.sims_fits[b], apers)
    sims_apflux= np.array(apy_table['aperture_sum'])

    sky_apflux= img_apflux - sims_apflux
    #sky_meas= t.obitractor[itrac].flux_z - t.simcat.zflux

    simcat_flux= t.simcat.get(b+'flux')
    trac_flux= t.obitractor[itrac].get('flux_'+b)
    trac_apflux=t.obitractor[itrac].get('apflux_'+b)[:,5]

    print('sky_apflux:',sky_apflux)
    #print('sky_meas:',sky_meas)
    print('sims_apflux:',sims_apflux)
    print('simcat_flux:',simcat_flux)
    print('obitrac_flux:',trac_flux)
    print('obitrac_apflux:',trac_apflux)
    print(sky_apflux + sims_apflux - trac_apflux)
    print(simcat_flux + sky_apflux - trac_apflux)
('band=', 'g')
('sky_apflux:', array([ 0.11485741,  0.03434163,  1.73282194,  0.73871268]))
('sims_apflux:', array([ 53.08364024,  53.13556137,  53.07870167,  53.12389901]))
('simcat_flux:', array([ 52.83479175,  52.83479175,  52.83479175,  52.83479175]))
('obitrac_flux:', array([ 76.45344543,  74.84441376,  78.35695648,  77.82241821], dtype=float32))
('obitrac_apflux:', array([ 53.10929108,  53.07685852,  54.82538605,  53.70322418], dtype=float32))
[ 0.08920657  0.09304447 -0.01386244  0.15938751]
[-0.15964192 -0.20772514 -0.25777236 -0.12971975]
('band=', 'r')
('sky_apflux:', array([ 0.11485741,  0.03434163,  1.73282194,  0.73871268]))
('sims_apflux:', array([ 53.08364024,  53.13556137,  53.07870167,  53.12389901]))
('simcat_flux:', array([ 52.83479175,  52.83479175,  52.83479175,  52.83479175]))
('obitrac_flux:', array([ 76.45344543,  74.84441376,  78.35695648,  77.82241821], dtype=float32))
('obitrac_apflux:', array([ 53.10929108,  53.07685852,  54.82538605,  53.70322418], dtype=float32))
[ 0.08920657  0.09304447 -0.01386244  0.15938751]
[-0.15964192 -0.20772514 -0.25777236 -0.12971975]
('band=', 'z')
('sky_apflux:', array([ 0.11485741,  0.03434163,  1.73282194,  0.73871268]))
('sims_apflux:', array([ 53.08364024,  53.13556137,  53.07870167,  53.12389901]))
('simcat_flux:', array([ 52.83479175,  52.83479175,  52.83479175,  52.83479175]))
('obitrac_flux:', array([ 76.45344543,  74.84441376,  78.35695648,  77.82241821], dtype=float32))
('obitrac_apflux:', array([ 53.10929108,  53.07685852,  54.82538605,  53.70322418], dtype=float32))
[ 0.08920657  0.09304447 -0.01386244  0.15938751]
[-0.15964192 -0.20772514 -0.25777236 -0.12971975]
In [137]:
a=np.ones((10,10))+34
a /= a.sum()
a.sum()
Out[137]:
0.99999999999999989
In [148]:
assert(len(itrac) == 4)
print(simcat[isim].zflux-obitractor[itrac].flux_z)
print(np.sqrt(1/obitractor[itrac].flux_ivar_z))
print(flux2mag(simcat[isim].zflux)-flux2mag(obitractor[itrac].flux_z))

[ 2.80877871  3.7907352   1.49074359  3.23818774]
[ 0.36161846  0.34976208  0.36395597  0.35320479]
[-0.08277248 -0.11326717 -0.04314349 -0.09600375]
In [102]:
# sizes?
print(simcat[isim].rhalf-obitractor[itrac].shapeexp_r)
[-0.04129598  0.0056425  -0.05493272 -0.01868245]
In [103]:
# Does extinction vary by this much? I doubt it...
mags= flux2mag(obitractor[itrac].flux_z)
print('measured mag variation=',mags - np.median(mags))
mags= flux2mag(obitractor[itrac].flux_z/obitractor[itrac].mw_transmission_z)
print('extinction variation=',mags - np.median(mags))
('measured mag variation=', array([-0.00476646,  0.02602196, -0.05817604,  0.00476646], dtype=float32))
('extinction variation=', array([-0.00471115,  0.02607346, -0.05801773,  0.00471115], dtype=float32))
In [104]:
np.std(obitractor.mw_transmission_z)
Out[104]:
6.3498112e-05

sim sources should have similar flux to real galaxy at center

In [105]:
# Diff b/w fake truth and measure real galaxy identical
print(simcat[isim].zflux-obitractor[ireal].flux_z)
print(flux2mag(simcat[isim].zflux)-flux2mag(obitractor[ireal].flux_z))
[-0.00122647 -0.00122647 -0.00122647 -0.00122647]
[  3.50646973e-05   3.50646973e-05   3.50646973e-05   3.50646973e-05]
In [106]:
# But diff b/w measure fake and measure real galaxy are large
print(obitractor[itrac].flux_z-obitractor[ireal].flux_z)
print(flux2mag(obitractor[itrac].flux_z)-flux2mag(obitractor[ireal].flux_z))
[-5.3830986  -6.30246735 -3.72509003 -5.67053986]
[ 0.164608    0.19539642  0.11119843  0.17414093]
In [110]:
obitractor.flux_z
Out[110]:
array([ 32.88297653,  38.26607513,  32.59553528,  34.54098511,
        31.96360779,   4.96016264,   2.12047195], dtype=float32)
In [107]:

In [113]:
plot3d(img_jpg[:,:,0])
In [114]:
plot3d(model_jpg[:,:,0])
In [22]:
fig,ax=plt.subplots()
plotImage().imshow(blobs,ax,qs=None)
In [27]:
ra=[174.24714700910047, 174.24715333044438, 174.24915898792247, 174.2451413528916]
dec=[24.326224120319477, 24.32988191327332, 24.328050866811363, 24.32805515572893]
np.mean(ra),np.mean(dec)
Out[27]:
(174.24715017008975, 24.328053014033273)
In [41]:
os.environ["LEGACY_SURVEY_DIR"]= IN_DIR

from legacypipe.survey import LegacySurveyData, wcs_for_brick
survey = LegacySurveyData()
brickinfo = survey.get_brick_by_name('1741p242')
brickwcs = wcs_for_brick(brickinfo)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-41-a294eb3f18c1> in <module>()
      4 survey = LegacySurveyData()
      5 brickinfo = survey.get_brick_by_name('1741p242')
----> 6 brickwcs = wcs_for_brick(brickinfo)

/home/kaylan/myrepo/legacypipe/py/legacypipe/survey.pyc in wcs_for_brick(b, W, H, pixscale)
    533     from astrometry.util.util import Tan
    534     pixscale = pixscale / 3600.
--> 535     return Tan(b.ra, b.dec, W/2.+0.5, H/2.+0.5,
    536                -pixscale, 0., 0., pixscale,
    537                float(W), float(H))

AttributeError: 'NoneType' object has no attribute 'ra'
In [28]:
brickwcs.radec2pixelxy(174.24715,24.32805)
Out[28]:
(True, 190.26994464884524, 2873.6879080454523)
In [38]:
_,x,y= brickwcs.radec2pixelxy(174.2471,24.3262)
x,y= int(x),int(y)+25
hw=100

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
plotImage().imshow(sliceImage(blobs,xslice=xslc,yslice=yslc),
                   ax,qs=None)
print(xslc,yslc)
(slice(90, 290, None), slice(2773, 2973, None))

DECam

In [27]:
# DR3: Correct Region?
from legacypipe.survey import LegacySurveyData, wcs_for_brick
dr= os.path.join(os.environ['HOME'],
                'myrepo/legacypipe/dr3_testcase')
os.environ["LEGACY_SURVEY_DIR"]= dr

survey = LegacySurveyData()
brick='2357p087'
brickinfo = survey.get_brick_by_name(brick)
brickwcs = wcs_for_brick(brickinfo)

ra,dec= 235.6799,8.8597
print(brickwcs.radec2pixelxy(ra,dec))

# Plot
_,x,y= brickwcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
img_jpg= readImage(os.path.join(dr,'legacysurvey-%s-image.jpg' % brick),
                   jpeg=True)
plotImage().imshow(sliceImage(img_jpg,xslice=xslc,yslice=yslc),
                   ax,qs=None)
print(x,y)
print(x-hw,x+hw,y-hw,y+hw)

for x,y in [(x-hw,y),(x+hw,y),(x,y+hw),(x,y-hw)]:
    print(brickwcs.pixelxy2radec(x,y))
(True, 2675.942971854415, 3307.9059576884274)
(2675, 3307)
(2575, 2775, 3207, 3407)
(235.68733510594618, 8.859635266642558)
(235.67260383081032, 8.85963274643333)
(235.67996820762602, 8.866911823308769)
(235.67997072695223, 8.852356330322968)
In [15]:
# TIM objects for each CCD
from legacypipe.decam import DecamImage
from legacypipe.survey import LegacySurveyData

ccds= fits_table(os.path.join(dr,'dr3_testcase_ccds.fits'))

kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=True, invvar=True)

survey = LegacySurveyData(ccds=ccds)
tims={}
#for col in ccds.get_columns():
#    if 'S' in str(ccds.get(col).dtype):
#        ccds.set(col,np.char.strip(ccds.get(col)))
for ccd in ccds:
    im = survey.get_image_object(ccd)
    tims[ccd.filter] = im.get_tractor_image(**kwargs)
Reading image slice: None
Reading image from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_g_v2/c4d_140818_235551_ooi_g_v2.fits.fz hdu 35
Reading inverse-variance from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_g_v2/c4d_140818_235551_oow_g_v2.fits.fz hdu 35
Reading data quality from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_g_v2/c4d_140818_235551_ood_g_v2.fits.fz hdu 35
Warning: un-setting SATUR for 3 pixels with SATUR and BADPIX set.
Merged spline sky model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky-merged/00349/decam-00349605.fits
Reading sky model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky/00349/00349605/decam-00349605-N5.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex-merged/00349/decam-00349605.fits
Reading PsfEx model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex/00349/00349605/decam-00349605-N5.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 2.09, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.114103
Galaxy norm 0.0902848326709
Reading image slice: None
Reading image from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_r_v2/c4d_140813_005917_ooi_r_v2.fits.fz hdu 35
Reading inverse-variance from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_r_v2/c4d_140813_005917_oow_r_v2.fits.fz hdu 35
Reading data quality from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_r_v2/c4d_140813_005917_ood_r_v2.fits.fz hdu 35
Warning: un-setting SATUR for 2 pixels with SATUR and BADPIX set.
Merged spline sky model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky-merged/00347/decam-00347323.fits
Reading sky model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky/00347/00347323/decam-00347323-N5.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex-merged/00347/decam-00347323.fits
Reading PsfEx model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex/00347/00347323/decam-00347323-N5.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.93, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.136238
Galaxy norm 0.102938388518
Reading image slice: None
Reading image from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_z_v2/c4d_140811_002713_ooi_z_v2.fits.fz hdu 35
Reading inverse-variance from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_z_v2/c4d_140811_002713_oow_z_v2.fits.fz hdu 35
Reading data quality from /home/kaylan/myrepo/legacypipe/dr3_testcase/images/decam/CP20140810_z_v2/c4d_140811_002713_ood_z_v2.fits.fz hdu 35
Warning: un-setting SATUR for 47 pixels with SATUR and BADPIX set.
Merged spline sky model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky-merged/00346/decam-00346606.fits
Reading sky model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/splinesky/00346/00346606/decam-00346606-N5.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex-merged/00346/decam-00346606.fits
Reading PsfEx model from /home/kaylan/myrepo/legacypipe/dr3_testcase/calib/decam/psfex/00346/00346606/decam-00346606-N5.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 2.42, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.106902
Galaxy norm 0.0885711598239
In [24]:
# x,y pixels for each ccd
hw=100
for b in ccds.filter:
    _,x,y= tims[b].wcs.wcs.radec2pixelxy(ra,dec)
    print('%s: ' % b,x-hw,x+hw,y-hw,y+hw)
('g: ', 318.70780228177273, 518.7078022817727, 1921.286605631049, 2121.286605631049)
('r: ', 367.16452653880515, 567.1645265388051, 1968.49883705385, 2168.49883705385)
('z: ', 357.70981507257636, 557.7098150725764, 1941.7682361872326, 2141.7682361872326)
In [43]:
# brick 0285m165, RA,Dec = (28.4194, -16.4362) in viewer looks like good 3 colors source
IN_DIR= os.path.join(os.environ['HOME'],
                     'myrepo/obiwan/tests/end_to_end',
                     'testcase_DR5_grz/'))
os.environ["LEGACY_SURVEY_DIR"]= IN_DIR

survey = LegacySurveyData()
brickinfo = survey.get_brick_by_name('0285m165')
brickwcs = wcs_for_brick(brickinfo)

ra,dec= 28.4194, -16.4362
brickwcs.radec2pixelxy(ra,dec)
Out[43]:
(True, 3177.4193025452664, 2676.786370753425)
In [67]:
brick='0285m165'
DATA_DIR= os.path.join(os.environ['HOME'],
                       'mydata',brick,'dr5')
dr5tractor= fits_table(os.path.join(DATA_DIR,
                      'tractor-%s.fits' % brick))
#ccds= fits_table(os.path.join(DATA_DIR,
#                 'dr5_3ccds.fits'))
#ccds= fits_table(os.path.join(DATA_DIR,
#                 'ccds_0285m165_dr5coadddir.fits'))
ccds= fits_table(os.path.join(IN_DIR,'survey-ccds-1.fits.gz'))

img_jpg= readImage(os.path.join(DATA_DIR,'legacysurvey-0285m165-image.jpg'),
                   jpeg=True)
In [46]:
_,x,y= brickwcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
plotImage().imshow(sliceImage(img_jpg,xslice=xslc,yslice=yslc),
                   ax,qs=None)
print(xslc,yslc)
(slice(3077, 3277, None), slice(2576, 2776, None))

Mosaic, 90Prime

In [4]:
# DR6: Correct Region?
from legacypipe.survey import LegacySurveyData, wcs_for_brick
dr= os.path.join(os.environ['HOME'],
                'myrepo/obiwan/tests/end_to_end/testcase_bassmzls_grz')
os.environ["LEGACY_SURVEY_DIR"]= dr

survey = LegacySurveyData()
brick='2176p330'
brickinfo = survey.get_brick_by_name(brick)
brickwcs = wcs_for_brick(brickinfo)

ra,dec= 217.5429, 33.0873
_,x,y= brickwcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100
print('brick box pixels: x=[%d,%d],y=[%d,%d]' % \
     (x-hw,x+hw,y-hw,y+hw))

for x,y in [(x-hw,y),(x+hw,y),(x,y+hw),(x,y-hw)]:
    print(brickwcs.pixelxy2radec(x,y))

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
img_jpg= readImage(os.path.join(dr,'legacysurvey-%s-image.jpg' % brick),
                   jpeg=True)
plotImage().imshow(sliceImage(img_jpg,xslice=xslc,yslice=yslc),
                   ax,qs=None)
---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-4-4d057a932720> in <module>()
      7 survey = LegacySurveyData()
      8 brick='2176p330'
----> 9 brickinfo = survey.get_brick_by_name(brick)
     10 brickwcs = wcs_for_brick(brickinfo)
     11

/srv/repos_for_docker/legacypipe/py/legacypipe/survey.py in get_brick_by_name(self, brickname)
   1161         Returns a brick (as one row in a table) by name (string).
   1162         '''
-> 1163         B = self.get_bricks_readonly()
   1164         I, = np.nonzero(np.array([n == brickname for n in B.brickname]))
   1165         if len(I) == 0:

/srv/repos_for_docker/legacypipe/py/legacypipe/survey.py in get_bricks_readonly(self)
   1140         '''
   1141         if self.bricks is None:
-> 1142             self.bricks = self.get_bricks()
   1143             # Assert that bricks are the sizes we think they are.
   1144             # ... except for the two poles, which are half-sized

/srv/repos_for_docker/legacypipe/py/legacypipe/survey.py in get_bricks(self)
   1133         uses a cached version.
   1134         '''
-> 1135         return fits_table(self.find_file('bricks'))
   1136
   1137     def get_bricks_readonly(self):

~/myinstall/astrometry/lib/python/astrometry/util/fits.py in fits_table(dataorfn, rows, hdunum, hdu, ext, header, columns, column_map, lower, mmap, normalize, use_fitsio, tabledata_class)
    665
    666         if fitsio:
--> 667             F = fitsio.FITS(dataorfn)
    668             data = F[hdunum]
    669             hdr = data.read_header()

/srv/py3_venv/lib/python3.5/site-packages/fitsio/fitslib.py in __init__(self, filename, mode, **keys)
    361                     create=1
    362
--> 363         self._FITS =  _fitsio_wrap.FITS(filename, self.intmode, create)
    364
    365     def close(self):

OSError: FITSIO status = 104: could not open the named file
failed to find or open the following file: (ffopen)
/root/myrepo/obiwan/tests/end_to_end/testcase_bassmzls_grz/survey-bricks.fits.gz

In [10]:
survey = LegacySurveyData(ccds=ccds)
im = survey.get_image_object(ccds[0])

kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=True, invvar=True)
im.get_tractor_image(**kwargs)
camera=-mosaic- camera.strip=-mosaic-
Reading image slice: None
Reading image from /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_ooi_zd_v1-CCD3.fits hdu 1
Reading weight map image /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_oow_zd_v1-CCD3.fits ext 1
Reading data quality image /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_ood_zd_v1-CCD3.fits ext 1
Remapping weight map for 00110480-CCD3
Merged spline sky model does not exist: /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/splinesky-merged/00110/mosaic-00110480.fits
Reading sky model from /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/splinesky/00110/00110480/mosaic-00110480-CCD3.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/psfex-merged/00110/mosaic-00110480.fits
Reading PsfEx model from /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/psfex/00110/00110480/mosaic-00110480-CCD3.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.81, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.145167
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-b9337c9d9c8d> in <module>()
      4 kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
      5               pixels=True, dq=True, invvar=True)
----> 6 im.get_tractor_image(**kwargs)

/home/kaylan/myrepo/legacypipe/py/legacypipe/image.pyc in get_tractor_image(self, slc, radecpoly, gaussPsf, pixPsf, hybridPsf, splinesky, nanomaggies, subsky, tiny, dq, invvar, pixels, constant_invvar)
    331         print('PSF norm:', tim.psfnorm)
    332         # Galaxy-detection norm
--> 333         tim.galnorm = self.galaxy_norm(tim)
    334         print('Galaxy norm', tim.galnorm)
    335         assert(tim.galnorm < tim.psfnorm)

/home/kaylan/myrepo/legacypipe/py/legacypipe/image.pyc in galaxy_norm(self, tim, x, y)
    437             y = h//2
    438         pos = tim.wcs.pixelToPosition(x, y)
--> 439         gal = SimpleGalaxy(pos, NanoMaggies(**{band:1.}))
    440         S = 32
    441         mm = ModelMask(int(x-S), int(y-S), 2*S+1, 2*S+1)

/home/kaylan/myrepo/legacypipe/py/legacypipe/survey.py in __init__(self, *args)
     91
     92     def __init__(self, *args):
---> 93         super(SimpleGalaxy, self).__init__(*args)
     94         self.shape = SimpleGalaxy.shape
     95

TypeError: super(type, obj): obj must be an instance or subtype of type
In [11]:
# TIM objects for each CCD
from legacypipe.mosaic import MosaicImage
from legacypipe.bok import BokImage
from legacypipe.survey import LegacySurveyData

ccds= fits_table(os.path.join(dr,'survey-ccds-1.fits.gz'))

kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=True, invvar=True)

survey = LegacySurveyData(ccds=ccds)
tims={}
#for col in ccds.get_columns():
#    if 'S' in str(ccds.get(col).dtype):
#        ccds.set(col,np.char.strip(ccds.get(col)))
for ccd in ccds:
    im = survey.get_image_object(ccd)
    tims[ccd.filter] = im.get_tractor_image(**kwargs)
camera=-mosaic- camera.strip=-mosaic-
Reading image slice: None
Reading image from /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_ooi_zd_v1-CCD3.fits hdu 1
Reading weight map image /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_oow_zd_v1-CCD3.fits ext 1
Reading data quality image /home/kaylan/mydata/testcase_dr6/test/images/mosaic/CP20170224/k4m_170225_115030_ood_zd_v1-CCD3.fits ext 1
Remapping weight map for 00110480-CCD3
Merged spline sky model does not exist: /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/splinesky-merged/00110/mosaic-00110480.fits
Reading sky model from /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/splinesky/00110/00110480/mosaic-00110480-CCD3.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/psfex-merged/00110/mosaic-00110480.fits
Reading PsfEx model from /home/kaylan/mydata/testcase_dr6/test/calib/mosaic/psfex/00110/00110480/mosaic-00110480-CCD3.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.81, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.145167
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-712d6a45f37e> in <module>()
     16 for ccd in ccds:
     17     im = survey.get_image_object(ccd)
---> 18     tims[ccd.filter] = im.get_tractor_image(**kwargs)

/home/kaylan/myrepo/legacypipe/py/legacypipe/image.pyc in get_tractor_image(self, slc, radecpoly, gaussPsf, pixPsf, hybridPsf, splinesky, nanomaggies, subsky, tiny, dq, invvar, pixels, constant_invvar)
    331         print('PSF norm:', tim.psfnorm)
    332         # Galaxy-detection norm
--> 333         tim.galnorm = self.galaxy_norm(tim)
    334         print('Galaxy norm', tim.galnorm)
    335         assert(tim.galnorm < tim.psfnorm)

/home/kaylan/myrepo/legacypipe/py/legacypipe/image.pyc in galaxy_norm(self, tim, x, y)
    437             y = h//2
    438         pos = tim.wcs.pixelToPosition(x, y)
--> 439         gal = SimpleGalaxy(pos, NanoMaggies(**{band:1.}))
    440         S = 32
    441         mm = ModelMask(int(x-S), int(y-S), 2*S+1, 2*S+1)

/home/kaylan/myrepo/legacypipe/py/legacypipe/survey.py in __init__(self, *args)
     91
     92     def __init__(self, *args):
---> 93         super(SimpleGalaxy, self).__init__(*args)
     94         self.shape = SimpleGalaxy.shape
     95

TypeError: super(type, obj): obj must be an instance or subtype of type
In [24]:
# x,y pixels for each ccd
hw=100
for b in ccds.filter:
    _,x,y= tims[b].wcs.wcs.radec2pixelxy(ra,dec)
    print('%s: ' % b,x-hw,x+hw,y-hw,y+hw)
('g: ', 318.70780228177273, 518.7078022817727, 1921.286605631049, 2121.286605631049)
('r: ', 367.16452653880515, 567.1645265388051, 1968.49883705385, 2168.49883705385)
('z: ', 357.70981507257636, 557.7098150725764, 1941.7682361872326, 2141.7682361872326)
In [43]:
# brick 0285m165, RA,Dec = (28.4194, -16.4362) in viewer looks like good 3 colors source
IN_DIR= os.path.join(os.environ['HOME'],
                     'myrepo/obiwan/tests/end_to_end',
                     'testcase_DR5_grz/'))
os.environ["LEGACY_SURVEY_DIR"]= IN_DIR

survey = LegacySurveyData()
brickinfo = survey.get_brick_by_name('0285m165')
brickwcs = wcs_for_brick(brickinfo)

ra,dec= 28.4194, -16.4362
brickwcs.radec2pixelxy(ra,dec)
Out[43]:
(True, 3177.4193025452664, 2676.786370753425)
In [67]:
brick='0285m165'
DATA_DIR= os.path.join(os.environ['HOME'],
                       'mydata',brick,'dr5')
dr5tractor= fits_table(os.path.join(DATA_DIR,
                      'tractor-%s.fits' % brick))
#ccds= fits_table(os.path.join(DATA_DIR,
#                 'dr5_3ccds.fits'))
#ccds= fits_table(os.path.join(DATA_DIR,
#                 'ccds_0285m165_dr5coadddir.fits'))
ccds= fits_table(os.path.join(IN_DIR,'survey-ccds-1.fits.gz'))

img_jpg= readImage(os.path.join(DATA_DIR,'legacysurvey-0285m165-image.jpg'),
                   jpeg=True)
In [46]:
_,x,y= brickwcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100

fig,ax=plt.subplots()
xslc=slice(x-hw,x+hw)
yslc=slice(y-hw,y+hw)
plotImage().imshow(sliceImage(img_jpg,xslice=xslc,yslice=yslc),
                   ax,qs=None)
print(xslc,yslc)
(slice(3077, 3277, None), slice(2576, 2776, None))

DECam

In [16]:
from legacypipe.decam import DecamImage
from legacypipe.survey import LegacySurveyData
In [80]:
kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=True, invvar=True)

survey = LegacySurveyData(ccds=ccds)
tims={}
#for col in ccds.get_columns():
#    if 'S' in str(ccds.get(col).dtype):
#        ccds.set(col,np.char.strip(ccds.get(col)))
for ccd in ccds:
    im = survey.get_image_object(ccd)
    tims[ccd.filter] = im.get_tractor_image(**kwargs)
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141206_043719_ooi_z_a1-S10.fits hdu 1
Reading inverse-variance from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141206_043719_oow_z_a1-S10.fits hdu 1
Reading data quality from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141206_043719_ood_z_a1-S10.fits hdu 1
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00384/decam-00384128.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00384/00384128/decam-00384128-S10.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00384/decam-00384128.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00384/00384128/decam-00384128-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.64, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.158392
Galaxy norm 0.113801322034
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141224_033245_ooi_r_a1-S10.fits hdu 1
Reading inverse-variance from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141224_033245_oow_r_a1-S10.fits hdu 1
Reading data quality from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141224_033245_ood_r_a1-S10.fits hdu 1
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00390/decam-00390924.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00390/00390924/decam-00390924-S10.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00390/decam-00390924.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00390/00390924/decam-00390924-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.49, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.178092
Galaxy norm 0.121385678047
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141223_020302_ooi_g_a1-S10.fits hdu 1
Reading inverse-variance from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141223_020302_oow_g_a1-S10.fits hdu 1
Reading data quality from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141223_020302_ood_g_a1-S10.fits hdu 1
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00390/decam-00390498.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00390/00390498/decam-00390498-S10.fits
Instantiating and subtracting sky model
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00390/decam-00390498.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00390/00390498/decam-00390498-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 2.05, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.123636
Galaxy norm 0.0966921909813
In [81]:
tim=tims['g']
tim.wcs.wcs.radec2pixelxy(ra,dec)
Out[81]:
(True, 918.0465384893796, 420.42404627704695)
In [82]:
_,x,y= tim.wcs.wcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100
xslc,yslc= slice(x-hw,x+hw),slice(y-hw,y+hw)

fig,ax=plt.subplots(1,3,figsize=(7,4))
for i,band in enumerate(['g','r','z']):
    plotImage().imshow(sliceImage(tims[band].getImage(),xslice=xslc,yslice=yslc).T,
                       ax[i])
In [87]:
def flux2mag(flux):
    return -2.5*np.log10(1e-9 * flux)

hw=30
xslc,yslc= slice(x-hw,x+hw),slice(y-hw,y+hw)

fig,ax=plt.subplots(1,3,figsize=(7,4))
for i,band in enumerate(['g','r','z']):
    img= sliceImage(tims[band].getImage(),xslice=xslc,yslice=yslc).T
    plotImage().imshow(img,ax[i])
    # Extract the nanomaggie flux for each band that we will add to the images
    print(band,'nanomags=',img.sum(),'mags=',flux2mag(img.sum()))
('g', 'nanomags=', 17.956371, 'mags=', 19.364453556756388)
('r', 'nanomags=', 32.756645, 'mags=', 18.711751458105404)
('z', 'nanomags=', 46.451378, 'mags=', 18.332503497939175)
In [72]:
# flux of this galaxy
from astrometry.libkd.spherematch import match_radec
_,itrac,d= match_radec(np.array([ra]), np.array([dec]), dr5tractor.ra, dr5tractor.dec,
                        5./3600.0,nearest=True)
itrac
Out[72]:
array([4962], dtype=int32)
In [74]:
trac=dr5tractor.copy()
trac.cut(itrac)
len(trac)
Out[74]:
1
In [76]:
trac.ra,trac.dec, ra,dec
Out[76]:
(array([ 28.4192355]), array([-16.43595489]), 28.4194, -16.4362)
In [88]:
# Tractor measurements on many imags should be close to above values for the three images
for band in ['g','r','z']:
    flux=trac.get('flux_'+band)
    print(band,'nanomags=',flux,'mags=',flux2mag(flux))
('g', 'nanomags=', array([ 16.06177902], dtype=float32), 'mags=', array([ 19.48551559], dtype=float32))
('r', 'nanomags=', array([ 31.8071003], dtype=float32), 'mags=', array([ 18.74369049], dtype=float32))
('z', 'nanomags=', array([ 52.83349609], dtype=float32), 'mags=', array([ 18.19272614], dtype=float32))
In [89]:
trac.ra,trac.dec,trac.type,flux2mag(trac.flux_g),flux2mag(trac.flux_r),flux2mag(trac.flux_z)
Out[89]:
(array([ 28.4192355]), array([-16.43595489]), array(['EXP '],
       dtype='|S4'), array([ 19.48551559], dtype=float32), array([ 18.74369049], dtype=float32), array([ 18.19272614], dtype=float32))
In [90]:
trac.shapeexp_r,trac.shapeexp_e1,trac.shapeexp_e2
Out[90]:
(array([ 2.4597683], dtype=float32),
 array([ 0.12245165], dtype=float32),
 array([ 0.23880503], dtype=float32))
In [104]:
# Lets circle ra,dec where we'll add galaxies to the r-band image
print(ra,dec)
tim= tims['r']
dd= 60 * 0.262 /3600 # pix*as/pix converted to deg
newra,newdec=[],[]
for dx,dy in zip([dd,-dd,0,0],[0,0,dd,-dd]):
    #newx,newy=x+dx,y+dy
    #ra,dec= tim.wcs.wcs.pixelxy2radec(newx,newy)
    #print(newx,newy,ra,dec)
    newra.append(ra + dx)
    newdec.append(dec + dy)
print(newra)
print(newdec)
_,newx,newy= tim.wcs.wcs.radec2pixelxy(newra,newdec)
print(newx)
print(newy)

_,x,y= tim.wcs.wcs.radec2pixelxy(ra,dec)
x,y= int(x),int(y)
hw=100
xslc,yslc= slice(x-hw,x+hw),slice(y-hw,y+hw)

fig,ax=plt.subplots(figsize=(3,3))
plotImage().imshow(sliceImage(tims['r'].getImage(),xslice=xslc,yslice=yslc),
                   ax)
# cirlces
plotImage().circles(newx,newy,ax,
                    xslice=xslc,yslice=yslc,
                    r_pixels=5./0.262,color='y')
print('IMAGE SLICE:',xslc,yslc)
(28.4194, -16.4362)
[28.423766666666666, 28.415033333333334, 28.4194, 28.4194]
[-16.4362, -16.4362, -16.431833333333334, -16.440566666666665]
[ 918.6980714   918.75985247  859.09916007  978.3570594 ]
[ 477.86830435  363.5315702   420.66734678  420.73262761]
('IMAGE SLICE:', slice(818, 1018, None), slice(320, 520, None))
In [125]:
# start with coadd/ccds
ccds= fits_table(os.path.join(DATA_DIR,
                 'ccds_0285m165_dr5coadddir.fits'))
In [111]:
kwargs = dict(pixPsf=True, splinesky=True, subsky=True, hybridPsf=True,
              pixels=True, dq=False, invvar=False)


survey = LegacySurveyData(ccds=ccds)
vals=[]
for band in ['g','r','z']:
    t= ccds[ccds.filter == band]
    for ccd in t:
        im = survey.get_image_object(ccd)
        tim = im.get_tractor_image(**kwargs)
        _,x,y= tim.wcs.wcs.radec2pixelxy(ra,dec)
        print(band,int(x),int(y))
        vals.append( (band,int(x),int(y)) )
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141223_020302_ooi_g_a1.fits.fz hdu 20
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00390/decam-00390498.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00390/00390498/decam-00390498-S10.fits
Instantiating and subtracting sky model
sig1 estimate: 0.00575166300675
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00390/decam-00390498.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00390/00390498/decam-00390498-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 2.10, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.123362
Galaxy norm 0.0972685759985
('g', 1814, 3205)
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141224_033245_ooi_r_a1.fits.fz hdu 20
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00390/decam-00390924.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00390/00390924/decam-00390924-S10.fits
Instantiating and subtracting sky model
sig1 estimate: 0.00801456716987
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00390/decam-00390924.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00390/00390924/decam-00390924-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.41, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.162455
Galaxy norm 0.120581835731
('r', 1783, 3215)
Reading image slice: None
Reading image from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/images/decam/NonDECaLS/CP20141215/c4d_141206_043719_ooi_z_a1.fits.fz hdu 20
Merged spline sky model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky-merged/00384/decam-00384128.fits
Reading sky model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/splinesky/00384/00384128/decam-00384128-S10.fits
Instantiating and subtracting sky model
sig1 estimate: 0.0346027164371
Merged PsfEx model does not exist: /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex-merged/00384/decam-00384128.fits
Reading PsfEx model from /home/kaylan/myrepo/obiwan/tests/end_to_end/testcase_DR5_grz/calib/decam/psfex/00384/00384128/decam-00384128-S10.fits
Using PSF model HybridPixelizedPSF: Gaussian sigma 1.53, Pix PixelizedPsfEx
-- creating constant PSF model for norms...
PSF norm: 0.137621
Galaxy norm 0.111473202692
('z', 1782, 3263)
[('g', 1814, 3205), ('r', 1783, 3215), ('z', 1782, 3263)]
In [114]:
print(vals)
dx=100
from collections import defaultdict
x0,x1,y0,y1= defaultdict(list),defaultdict(list),defaultdict(list),defaultdict(list)
for b,xc,yc in vals:
    x0[b] += [=xc-100,xc+100,yc-100,yc+100)
[('g', 1814, 3205), ('r', 1783, 3215), ('z', 1782, 3263)]
('g', 1714, 1914, 3105, 3305)
('r', 1683, 1883, 3115, 3315)
('z', 1682, 1882, 3163, 3363)
In [122]:
IN_DIR= os.path.join(os.environ['HOME'],
                     'myrepo/obiwan/tests/end_to_end',
                     'testcase_DR5_grz_200x200/')
g=fitsio.FITS(os.path.join(IN_DIR,'images/decam/NonDECaLS/CP20141215/',
                             'c4d_141223_020302_ooi_g_a1-S10.fits'))[1].read()
r=fitsio.FITS(os.path.join(IN_DIR,'images/decam/NonDECaLS/CP20141215/',
                             'c4d_141224_033245_ooi_r_a1-S10.fits'))[1].read()
z=fitsio.FITS(os.path.join(IN_DIR,'images/decam/NonDECaLS/CP20141215/',
                             'c4d_141206_043719_ooi_z_a1-S10.fits'))[1].read()
g.shape
Out[122]:
(200, 200)
In [123]:
fig,ax=plt.subplots(1,3,figsize=(6,3))
plotImage().imshow(g,ax[0])
plotImage().imshow(r,ax[1])
plotImage().imshow(z,ax[2])

It worked!

6’’ separation + testcase_DR_z

In [5]:
OUT_DIR= os.path.join(os.environ['HOME'],
                       'myrepo/obiwan/tests/end_to_end',
                       'out_testcase_DR5_z/elg/174/1741p242/rs0/')
IN_DIR= os.path.join(os.environ['HOME'],
                       'myrepo/obiwan/tests/end_to_end',
                       'testcase_DR5_z/')

simcat= fits_table(os.path.join(OUT_DIR,'obiwan/simcat-elg-1741p242.fits'))
obitractor= fits_table(os.path.join(OUT_DIR,'tractor/tractor-1741p242.fits'))

blobs= fitsio.FITS(os.path.join(OUT_DIR,'metrics/blobs-1741p242.fits.gz'))[0].read()

img_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-image.jpg'),
                   jpeg=True)
model_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-model.jpg'),
                     jpeg=True)
resid_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-resid.jpg'),
                     jpeg=True)
In [6]:
len(simcat),len(obitractor),blobs.shape
Out[6]:
(4, 3, (200, 200))

allblobs == False

Bug in Tractor, 3/5 sources modeled because they are apparently too close for the algorithm!

In [7]:
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(blobs,ax[0,0],qs=None)
plotImage().imshow(img_jpg,ax[0,1],qs=None)
plotImage().imshow(model_jpg,ax[1,0],qs=None)
plotImage().imshow(resid_jpg,ax[1,1],qs=None)

allblobs == True

Same Bug in Tractor, even when allblobs (there is only 1 blob anyway)

In [8]:
OUT_DIR= os.path.join(os.environ['HOME'],
                       'myrepo/obiwan/tests/end_to_end',
                       'out_testcase_DR5_z_allblobs/elg/174/1741p242/rs0/')
IN_DIR= os.path.join(os.environ['HOME'],
                       'myrepo/obiwan/tests/end_to_end',
                       'testcase_DR5_z_allblobs/')

simcat= fits_table(os.path.join(OUT_DIR,'obiwan/simcat-elg-1741p242.fits'))
obitractor= fits_table(os.path.join(OUT_DIR,'tractor/tractor-1741p242.fits'))

blobs= fitsio.FITS(os.path.join(OUT_DIR,'metrics/blobs-1741p242.fits.gz'))[0].read()

img_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-image.jpg'),
                   jpeg=True)
model_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-model.jpg'),
                     jpeg=True)
resid_jpg= readImage(os.path.join(OUT_DIR,'coadd/legacysurvey-1741p242-resid.jpg'),
                     jpeg=True)
In [9]:
len(simcat),len(obitractor),blobs.shape
Out[9]:
(4, 3, (200, 200))
In [10]:
fig,ax=plt.subplots(2,2,figsize=(6,6))
plotImage().imshow(blobs,ax[0,0],qs=None)
plotImage().imshow(img_jpg,ax[0,1],qs=None)
plotImage().imshow(model_jpg,ax[1,0],qs=None)
plotImage().imshow(resid_jpg,ax[1,1],qs=None)